function y = pk(m, rho, B, k, p0)

if ( 0 <= k && k <= m )
    y = p0 * ( ((m * rho)^k) / factorial(k));
else
    y = p0 * ( ((m * rho)^m) / factorial(m)) * ((rho)^(k-m));
end

end